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Abstract 

We present all-particle primary cosmic-ray energy spectrum in the 3 • 10 6 — 2 • 
10 8 GeV energy range obtained by a multi-parametric event-by-event evaluation 
of the primary energy. The results are obtained on the basis of an expanded EAS 
data set detected at mountain level (700 g/cm 2 ) by the GAMMA experiment. The 
energy evaluation method has been developed using the EAS simulation with the 
SIBYLL interaction model taking into account the response of GAMMA detectors 
and reconstruction uncertainties of EAS parameters. Nearly unbiased (< 5%) energy 
estimations regardless of a primary nuclear mass with an accuracy of about 15 — 10% 
in the 3 • 10 6 — 2 • 10 8 GeV energy range respectively are attained. 

An irregularity ('bump') in the spectrum is observed at primary energies of 
~ 7.4 • 10 7 GeV. This bump exceeds a smooth power-law fit to the data by about 
4 standard deviations. Not rejecting stochastic nature of the bump completely, we 
examined the systematic uncertainties of our methods and conclude that they can- 
not be responsible for the observed feature. 
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1 Introduction 



Study of the fine structure in the primary energy spectrum is one of the 
most important tasks in the very high energy cosmic ray experiments pQ . Com- 
monly accepted values of the all-particle energy spectrum indexes of —2.7 and 
—3.1 before and after the knee are an average and may not reflect the real 
behavior of the spectrum particularly after the knee. It is necessary to pay spe- 
cial attention to the energy region of (1 — 10) ■ 10 7 GeV, where experimental 
results have been very limited up to now. Irregularities of the energy spectrum 
in this region were observed a long time ago. They can be seen from energy 
spectrum obtained more than 20 years ago with AKENO experiment [2] as 
well as in later works of the GAMMA p] and TUNKA [3] experiments. At 
the same time the large statistical errors did not allow to discuss the reasons 
of these irregularities. 

On the other hand results of many experiments on the study of EAS charge 
particle spectra, the behavior of the age parameter and muon component char- 
acteristics point out that the primary mass composition at energies above the 
knee becomes significantly heavier. Based on these indications, additional in- 
vestigations of the fine structure of the primary energy spectrum at (1 — 10) • 10 7 
GeV have an obvious interest. 

There are two ways to obtain the primary energy spectra using detected 
extensive air showers (EAS). The first way is a statistical method, which un- 
folds the primary energy spectra from the corresponding integral equation set 
based on a detected EAS data set and the model of the EAS development in 
the atmosphere [51lo|TfSf5] . The second method is based on an event-by-event 
evaluation [2llO|lll|ll2|[T3"] of the primary energy of the detected EAS with pa- 
rameters q = q(N e , N^, Nh, s,9) using parametric E = /(q) [2|10|I11|I13] or 
non-parametric [12J energy estimator previously determined on the basis of 
shower simulations in the framework of a given model of EAS development. 

Here, applying a new event-by-event parametric energy evaluation E = 
/(q), the all-particle energy spectrum in the knee region is obtained on the 
basis of the data set obtained with the GAMMA EAS array [?1l8ll9lfTT] and a 
simulated EAS database obtained using the SIBYLL [T3] interaction model. 
Preliminary results have been presented in [TOlTT] . 



2 GAMMA experiment 

The GAMMA installation |7,8. 9|TT] is a ground based array of 33 surface 
detection stations and 150 underground muon detectors, located on the south 
side of Mount Aragats in Armenia. The elevation of the GAMMA facility is 
3200 m above sea level, which corresponds to 700 g/cm 2 of atmospheric depth. 
A diagrammatic layout of the array is shown in Fig. 1. 
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Fig. 1. Diagrammatic layout of the GAMMA facility. 



The surface stations of the EAS array are arranged in 5 concentric circles 
of ~20, 28, 50, 70 and 100 m radii, and each station contains 3 plastic scintil- 
lation detectors with the dimensions of 1x1x0. 05 m 3 . Each of the central 
9 stations contains an additional (the 4th) small scintillator with dimensions 
of 0.3 x 0.3 x 0.05 m 3 (Fig. 1) for high particle density 10 2 particles/m 2 ) 
measurements. 

A photomultiplier tube is placed on the top of the aluminum casing cover- 
ing each scintillator. One of three detectors of each station is viewed by two 
photomultipliers, one of which is designed for fast timing measurements. 

150 underground muon detectors ('muon carpet') are compactly arranged 
in the underground hall under 2.3 kg/cm 2 of concrete and rock. The scintil- 
lator dimensions, casings and photomultipliers are the same as in the EAS 
surface detectors. 

The shower size thresholds of the 100% shower detection efficiency are equal 
to N ch = 3 • 10 5 and N ch = 5 • 10 5 at the EAS core location within R < 25 m 
and R < 50 m respectively 7. 8HTT1 . 

The time delay is estimated by the pair-delay method [15J to give the time 
resolution of about 4 — 5 ns. The EAS detection efficiency (P<j) and corre- 
sponding shower parameter reconstruction errors are equal to: = 100%, 
A6 ~ 1.5°, AN ch /N ch ~ 0.05 - 0.15, As ~ 0.05, Ax and Ay ~ 0.7 - 1 m. 
The reconstruction errors of the truncated muon shower sizes for P M < 50 m 
from the shower core are equal to AN^/N^ ~ 0.2 — 0.35 at ~ 10 5 — 10 3 
respectively 8 ! )|TTj . 
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3 Event-by-event energy estimation 



3.1 Key assumptions 

Suppose that E\ = /(q) is an estimator of energy Eq of unknown primary 
nuclei which induced showers with the detected parameter q = q(N ch , N^, s, 6). 
Then the expected all-particle energy spectrum F(E{) is defined by 



where ^s(Eq) are the energy spectrum of primary nuclei and W(Eq, E\) are 
the corresponding (Eq, E{) transformation probability density function. 

If $s(Eq) oc Eq 1 and W(E ,Ei) are the log-normal distributions with 
S = Ei/Eq and a parameters, the expression (1) has the analytic solution 
for the expected spectrum of the energy estimator [i~6] : 



It is seen that evaluation of energy spectrum ^s(Eq) from (2) is possible to 
perform only at a priori known 7, 5 and a parameters and spectral slope 
(7) of detected energy spectra F(E\] coincides with spectral slope of primary 
energy spectra ^s(Eq). The values of 5 and a may depend on the primary 
energy (Eq) and mass of primary nuclei (A) from which the all-particle energy 
spectrum $s(E ) = J2a ^a(Eq) is consisted of. In this case, the expression (1) 
is unfolded numerically and the slope of detected energy spectrum can differ 
from primary energy spectrum. 

For example, the dependence o~(Eq) = aln (Eq/Eq^^) +b at \a\ < 0.1 leads 
to the numerical solutions which can be approximated by the expression (2) 
replacing a with o~(Eq) — a^/7. The corresponding approximation errors is 
about 2 — 5% in the energy range of E min — SOO-Emm and 7 ~ 2.3 — 3.2. 

However, the evaluation of energy spectra can be simplified provided 




(1) 




(2) 



j(Eq) ~ 7 ± A7, 



(3) 



5(E )~5 A (E 1 )~5 = 1±A6(E 1 ), 



(4) 



<t(Eq) 



a±Aa 



(5) 
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are satisfied for given energy range of E\. Then, the all-particle energy spec- 
trum can be evaluated from 

3(£ ) = F{E X ) exp ( - ((7 ~ 1V)2 ) . (6) 



The corresponding error of evaluation (6) with approaches (3-5) is determined 
by a sum of the statistical errors AF(Ei) and systematic errors rj due to 
approaches being used: 

(f) 2 ^) 2 + ^ 

where the systematic relative errors rj is 

r? 2 * (A*( 7 - l)) 2 + [((7(7 - 1)) 2 (^ + ^L)] 2 . (7) 



The values of 7, 5 and a parameters from expressions (3-7) and correspond- 
ing uncertainties A7, A5 and Act essential for the reconstruction of primary 
energy spectrum using the GAMMA facility EAS data and approach (6) are 
considered in Sections 3.2-3.4 below. 



3.2 Uncertainty of spectral slope 



The results of different experiments [9|17|ll8|ll9j and theoretical predictions 
[2"U|21f22j indicate that the elemental energy spectra can be presented in the 
power law form 



where 7 = 71 ~ 2.7ljj;f at E < E k (A) and 7 = 72 ~ 3.15±8:os at E > E k {A). 
It is also accepted that the mass spectra of primary nuclei can be divided 
into separate nuclear groups and below, as in [9], just 4 nuclear species (H, 
He, O-like and Fe-like) are considered. Dependence of the knee energy E k on 
the primary nuclei type assumed to be either rigidity dependent, E k — ZEz=\ 
[9l20f21|l22j or A-dependent [9"f2"3"] . E k = AEa=i, where Z and A are the charge 
and mass of primary nuclei correspondingly. 

As a result, the all-particle energy spectrum J2 a ^a{E) slowly changes its 
slope and can be roughly approximated by a power law spectrum with power 
index 7 ~ 2.7 at E < 3 • 10 6 GeV, 7 ~ 2.9 at 3 • 10 6 < E < 10 7 GeV and 
7 ~ 3.1 at E > 10 7 GeV. This presentation of the all-particle spectrum agrees 
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with world data [23J in the A7 ~ 0.1 range of uncertainty and energy interval 
10 6 < E < 2 ■ 10 8 GeV. 

The values of A5a(E) and cta(-E'o) parameters are presented in Section 3.4 
and depend on efficiency of energy estimator E x = f(N ch , N^, s, 9). 

Notice that it follows from the expression (7) that for o ~ 0.1 — 0.15 and 
Act = 0.03 the contribution of A7 in the systematic errors (7) is negligible and 
the difference of all-particle spectra evaluated by expression (6) for 7 = 2.7 
and 7 = 3.1 is less than 2% at ct ~ 0.15. 



3.3 The simulated EAS database 

To obtain the parametric representation for unbiased (8 ~ 1) energy esti- 
mator Ei of the primary energy Eq we simulated showers database using the 
CORSIKA(NKG) EAS simulation code [23] with the SIBYLL [H] interaction 
model for H, He, O and Fe primary nuclei. 

Preliminary, the showers simulated with NKG mode of CORSIKA code 
for each of the primary nuclei were compared with the corresponding simu- 
lations using EGS mode of CORSIKA [23] taking into account the detector 
response, contribution of EAS 7-quanta and shower parameter reconstruc- 
tion uncertainties. Simulated statistics were equal to 200 events for each of 
primary nuclei with log-uniform primary energy distribution in the range of 
2 • 10 6 — 10 8 GeV. Using the threshold energy of shower electrons (positrons) 
for NKG mode at observation level as a free parameter (the same as it was 
performed in [9]), the biases 5(N ch ,A) = (N ch (NKG) / N ch (EGS) ) - 1) and 
S S (A) = s(NKG) — s(EGS) were minimized for all simulated primary nuclei 
(if, He,Q,Fe). 

Applied method of calibration of the NKG mode of CORSIKA for the 
GAMMA EAS array differed from [9J only by the expanded range of selected 
shower core coordinates (R < 50m) and zenith angles 9 < 45°. The obtained 
biases of shower size S(N c h) and age parameter s in the range of statistical 
errors (< 5%) agreed with data [9]. The values of S(N ch ) were used further for 
correction of the shower size obtained by NKG simulation mode. 

The simulated primary energies (E ) for shower database were distributed 
according to a power law spectrum I(Eq) oc Eq 1 ' 5 with J\f = 2 ■ 10 4 total 
number of detected (N c h > 5 ■ 10 5 , R < 50m) and reconstructed showers for 
each primary nucleus. The energy thresholds of primary nuclei were set as 
Eo, m i n (A) = 10 6 GeV and -E max = 5 ■ 10 8 GeV. The simulated showers had 
core coordinates distributed uniformly within the radius of R < 75 m and 
zenith angles 9 < 45°. 

The reconstruction errors CT(lnA c / l ) of shower size N c h are presented in 
Fig. 2 for different primary nuclei and different zenith angles. The right and 
left ends of diagonals of the rectangular in Fig. 2 show the average primary 
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EAS size, N„ 



Fig. 2. Shower size reconstruction errors for different primary nuclei (p, He, O, Fe) 
and zenith angles (6 < 45° and 30° < 9 < 45°). The right and left ends of diagonals 
of the rectangular show the average primary energies (Eq) and corresponding shower 
sizes computed for the primary proton and Iron nuclei respectively. 



energies (in units of GeV) responsible for corresponding shower sizes for the 
primary proton and Iron nuclei respectively. 

All EAS muons with energies of > 4 GeV at GAMMA observation level 
have passed through the 2.3 kg/cm 2 of rock to the muon scintillation carpet 
(the underground muon hall, Fig. 1). The muon ionization losses and electron 
(positron) accompaniment due to muon electromagnetic and photonuclear in- 
teractions in the rock are taken into account using the approximation for equi- 
librium accompanying charged particles obtained from preliminary simulations 
with the FLUKA code [26J in the 0.005 — 20 TeV muon energy range. The 
resulting charged particle accompaniment per EAS muon in the underground 
hall is equal to 0.06 ±0.01 (100%e) and 11.0 ±1.5 (98.5%e, 1.4%h, 0.04%//) at 
muon energies 0.01 TeV and 10 TeV respectively. 

Due to absence of saturation in the muon scintillation carpet, the recon- 
struction errors (AlnA^) of truncated muon size are continuously de- 
creasing with increasing muon truncated sizes in the range 10 3 < < 10 5 . 
Corresponding magnitudes of reconstruction errors for primary protons and 
Iron nuclei were equal to A(lnA^ iP ) ~ 0.35,0.18,0.15 and A(lnA^ i i? e ) ~ 
0.38, 0.19, 0.10 for EAS muon truncated sizes ~ 10 3 , 10 4 , 10 5 respectively. 

Fluctuations of the shower size for given primary energies E QA = 10 6 , 10 7 , 10 8 
GeV and cos6> < 0.95 were equal to a A = p (N ch ,E ) ~ 0.20,0.14,0.10 and 
a AsFe {N ch: E ) ~ 0.16,0.13,0.08 respectively. 

Corresponding fluctuations of muon truncated size were equal to a A =p(N fJ/ , E ) 
~ 0.25,0.23,0.2 and a A =Fe(N^, E ) ~ 0.13,0.10,0.08. For zenith angles of pri- 
mary nuclei 45° > 9 > 30°, the fluctuations are increased about 1.5 — 2 times 
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Table 1 



Correlation coefficients p(q, InE'o) and /o(q, In A) between shower parameter q = 
q(N c h, N^, s) and primary energy (In E$) and nuclei mass In A for two zenith angular 
intervals. 



q 


lnE ,(9 < 10°) 


]nEo,(9 < 45°) 


In A, (9 < 10°) 


In A, (9 <45°) 


In N ch 
lniV^ 

s 


0.986 ± 0.001 
0.978 ± 0.001 
-0.029 ±0.013 


0.954 ± 0.0004 
0.969 ± 0.0003 
-0.02 ±0.004 


0.013 ±0.013 
0.139 ±0.012 
0.018 ±0.013 


0.007 ± 0.004 
0.132 ±0.004 
0.015 ± 0.004 



due to the aging of detected showers. 

The 4 x 2-10 4 EAS simulated events with reconstructed N C h, N^R < 50m), 
s and 9 shower parameters for the E and A parameters of primary nuclei made 
up the simulated EAS database. 



3.4 Energy estimator 



The event-by-event reconstruction of primary all-particle energy spectrum 
using the GAMMA facility is mainly based on high correlation of primary 
energy E and shower size (N ch ). The shower age parameter (s) zenith angle 
(9) and muon truncated shower size (A^) have to decrease the unavoidable 
biases of energy evaluations due to abundance of different primary nuclei. In 
Table 1 the correlation coefficients p(q, lnE ) and p(q, In A) between shower 
parameters N^, N^, s and primary energy (E ) and mass of primary nuclei 
(A = 1, 4, 16, 56) are presented. 

Parametric representation for the energy estimator In E 1 ~ f(a\N ch , iV M , s, 9) 
we obtained by minimizing % 2 

X 2 = 1211 ( lnE °' Ai - lnEl ^ 2 (9) 

A i=l °" 2 



with respect to a = a(ai, 02, . . . , a p ) for different empirical functions / (a) with 
a different number (p) of unknown parameters. The values of A, Eq and cor- 
responding reconstructed shower parameters N c h, N^, s and 9 for estimation 
of Ei were taken from simulated EAS database (Section 3.3). 

The best energy estimations as a result of the minimization (9) were achieved 
for the 7-parametric (j> = 7) fit: 

a 2 \/i CZ5 

In Ei = aix H h a 3 + a^c + -. + a-jye , (10) 

c [x — a^y) 
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Table 2 

Approximation parameters a±, . . . ,07 of primary energy evaluation (10) obtained 
from ^-minimization (9) for the SIBYLL interaction model, a = 0.14 and 
X 2 mi J n d.f. ^ 1- 



ai 


0.2 


03 


04 


a 5 


a 6 


a 7 


1.030 


3.641 


-5.743 


2.113 


6.444 


1.200 


-0.045 




6 7 8 

Log( E [ GeV ] ) 



Fig. 3. Mean biases versus energy E = Eq and E = E\ for the primary proton (p) 
and iron (Fe) nucleus and the uniformly mixed p, He,0, Fe compositions (All). 



where x = hiN c h, y = InN^R < 50m), c = cos 8, s is the shower age and 
energy E\ is in GeV. The values of a±, . . . , aj parameters are shown in Table 2 
and were derived at a = 0.14 and Xmm/ n d.f. — 1> where the number of degrees 
of freedom Ud.f. = 8 • 10 4 . The expected errors Aa±, . . . , Aaj of corresponding 
parameters were negligibly small (< 5%) due to very high values of nd.f.. 

The corresponding average biases versus energies (E = E and E = E{) for 
the primary proton (p), iron (-Fe) nucleus and uniformly mixed p, He,0, Fe 
composition are presented in Fig. 3 (symbols). The boundary lines correspond 
to approximations A5 ~ b/ E/10 6 GeV , where b ~ 0.10 and b ~ —0.17 for the 
upper and lower limits respectively. The shaded area corresponds to b ~ 0.09 
and b ~ —0.15 and were used to estimate the errors according to (7) for the 
reconstruction of all-particle energy spectrum (Section 4). 

The dependence of standard deviations a(E ) of systematic errors of energy 
evaluations (10) on primary energy E is presented in Fig. 4 for 4 primary 
nuclei and uniformly mixed composition with equal fractions of p, He, O and 
Fe nuclei. The results for the uniformly mixed composition with the shower 
core selection of R < 25 m [10J are presented in Fig. 4, as well. It is seen that 
the value of a = 0.14 responsible for x 2 — 1 (expression (9)) with uncertainty 
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Fig. 4. Errors of the energy estimator (6) versus primary energy Eq for 4 primary 
nuclei and uniformly mixed (Alt) composition. The empty rhombic symbols are 
taken from our previous data [TO] computed for the mixed composition and shower 
core selection criteria R < 25m. 




Fig. 5. Eq — Ei (in units of GeV) scatter plots for four (p, He, O, Fe) primary nuclei. 
The white lines show the corresponding Eq = E\ dependence. 
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Fig. 6. Distribution functions of errors for different primary nuclei (symbols) and 
uniformly mixed composition (symbols and solid line). 



Act ~ 0.03 (expression (5)) encloses the (Ja(Eo) data presented in Fig. 4. 

Such high accuracies of the energy evaluation regardless of primary nuclei 
is a consequence of the high mountain location of the GAMMA facility (700 
g/cm 2 ), where the correlation of primary energy with detected EAS size is 
about 0.95 - 0.99 (Table 1). 

The E — Ei scatter plot of simulated primary energy E Q and estimated 
energy Ei(N ch , N^, s, 9) according to expression (10) and Table 2 are shown in 
Fig. 5. The corresponding distributions of energy errors or the kernel function 
Wa(Eo, E\\5a, &a) of integral equation (1) for different primary nuclei and 
uniformly mixed composition are presented in Fig. 6 (symbols). The average 
values 5a(E ) and standard deviations <ja(E ) of these distributions depending 
on energy of primary nucleus (A) are presented in Figs. 3,4. The dashed line is 
an example of log-normal distribution with 5 and a parameters corresponding 
to the uniformly mixed composition. 

It is seen, that the errors can be described by the log-normal distributions 
and the key assumptions (3-5) are approximately valid. 

The test of applied approaches (expressions (3-6), Section 3.2) for the re- 
construction of all-particle primary energy spectrum was carried out by the 
direct folding of the power law energy spectrum ^s(E Q ) = dE /dE (expression 
(8)) with the log-normal kernel function W(E , Ei\8(E ),a(E )) according to 
expression (1) for primary proton and Iron nucleus. The values of S(E ) and 
<j(Eo) were derived from the log-parabolic interpolations of corresponding de- 
pendencies presented in Figs. 3,4. The event-by-event reconstructed energy 
spectrum dFi/dE\ (the left hand side of expression (1)) was obtained from 
expression (6) using approaches (3-5) with a = 0.14 and Act = 0.03. The 
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Fig. 7. Discrepancies of initial (dFo/dEo) and reconstructed (dF\/dEi) energy 
spectra (symbols) for the different primary nuclei and spectral indices of initial 
spectrum. The shaded area shows the expected errors according to expression (7). 
The star symbols are the spectral discrepancies for a pulsar component (Section 5). 



boundary lines of shaded area in Fig. 3 were used as estimations of uncertain- 
ties A<J(f?i) of condition (4). In Fig. 7 the values of (dF 1 /dE 1 ) / (dF /dE ) 
are presented (symbols) for primary proton and Iron nuclei and different "un- 
known" spectral indices of primary energy spectra (8) with the rigidity de- 
pendent knee at — 3 • 10 6 GeV. The shades area is the expected errors 
computed according to expression (7). It is seen, that all spectral discrep- 
ancies are practically covered by the expected errors according to expression 
(7). The star symbols in Fig. 7 represent the discrepancies of singular spectra 
with knee at energy 7.4 • 10 7 GeV described in Section 5. 



4 All-particle primary energy spectrum 



The EAS data set analyzed in this paper has been obtained for 5.63 • 10 sec 
of live run time of the GAMMA facility, from 2004 to 2006. Showers to be ana- 
lyzed were selected with the following criteria: N ch > 5-10 5 , R < 50 m, 9 < 45°, 
0.3 < s < 1.6, ^{N ch )/m < 3 and x 2 (V r M )/m < 3 (where m is the number of 
scintillators with non-zero signal), yielding a total data set of ~ 7- 10 5 selected 
showers. The selected measurement range provided the 100% EAS detection 
efficiency and similar conditions for the reconstruction of showers produced 
by primary nuclei H, He, . . . , Fe with energies 3 • 10 6 < E < (2 — 3) • 10 s 
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cos© 

Fig. 8. Detected zenith angular distributions for different energy thresholds (sym- 
bols). The lines are corresponding simulated isotropic distributions with the same 
statistics. 



GeV. The upper energy limit is determined from Fig. 4, where the saturation 
of surface scintillators in the shower core region begins to be significant. 

The independent test of energy estimates can be done by the detected 
zenith angle distributions which have to be isotropic for different energy 
thresholds. In Fig. 8 the corresponding detected distributions (symbols) are 
compared with statistically equivalent simulated isotropic distributions (lines). 
The agreement of detected and simulated distributions at E > 3 ■ 10 6 GeV 
gives an additional support to the consistency of energy estimates in the 
whole measurement range. The anisotropic spectral behavior at low energies 
(E ~ (1 — 3) • 10 6 GeV) is explained by the lack of heavy nuclei at larger zenith 
angles in the detected flux due to the applied shower selection criteria. 

Using the aforementioned unbiased (< 5%) event-by-event method of pri- 
mary energy evaluation (10), we obtained the all-particle energy spectrum. 
Results are presented in Fig. 9 (filled circle symbols, GAMMA07) in compar- 
ison with the same spectra obtained by the EAS inverse approach (line with 
shaded area, GAMMA06) from PI9"] and our preliminary results (point-circle 
symbols, GAMM A05) obtained using the 7-parametric event-by-event method 
with the shower core selection criteria R < 25m and 6» < 30° [10]. 

It follows from our preliminary data [TOirTT] . that the all-particle energy 
spectrum derived by event-by-event analysis with the multi-parametric en- 
ergy estimator (Section 3) depends only slightly on the interaction model 
( QGSJET01 [2H] or SIBYLL2.1 [H] ) and thereby, the errors of obtained 
spectra are mainly determined by the sum of statistical and systematic errors 
(7) presented in Fig. 9 by the dark shaded 

Shower size detection threshold effects distort the all-particle spectrum in 
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Fig. 9. All-particle energy spectrum in comparison with the results of EAS inverse 
approach |6|9j and our preliminary data TO] . The AKENO, Tibet-Ill, Fly's Eye 
Stereo, Hires/MIA and Hires-2 data were taken from [2127128 29 30J respectively. 



the range of E < (2 — 2.5) • 10 6 GeV depending on the interaction model and 
determine the lower limit -Emm = 3 • 10 6 GeV of the energy spectrum in Fig. 9 
whereas the upper limit of the spectrum -E max ~ (2 — 3) • 10 8 GeV is deter- 
mined by the saturation of our shower detectors which begins to be significant 
at E p > 2- 10 8 GeV and E Fe > 4- 10 s GeV (see Fig. 4) for primary proton and 
Fe nuclei. The range of minimal systematic errors and biases is (1 — 10) ■ 10 7 
GeV, where about 13% and 10% errors were attained (Fig. 3,4) for primary 
H and Fe nuclei respectively. 

In Table 3 the numerical values of the obtained all-particle energy spectrum 
are presented along with statistical, total upper and lower errors according to 
(7) and corresponding number of detected events. The energy spectra for low 
energy region (the first four lines) were taken from our data [10] for the EAS 
selection criteria R < 25m and 9 < 30°. 

The obtained energy spectrum agrees within errors with the KASCADE 
[6], AKENO [2] and Tibet-Ill [27J data both in the slope and in the abso- 
lute intensity practically in the whole measurement range. Looking at the 
experimental points we can unambiguously point out at the existence of an 
irregularity in the spectrum at the energy of (6 — 8) • 10 7 GeV. As it is seen 
from Figs 3 and 4, the energy estimator (10) has minimal biases (~ 4 — 5%) 
and errors (~ 0.09 — 0.12) at this energy. With these errors the obtained bump 
has an apparently real nature. If we fit all our other points in the (5 — 200) • 10 6 
GeV energy range by a smooth power-law spectrum, the bin at 7.4 • 10 7 GeV 
exceeds this smooth spectrum by 4.0 standard deviations. The exact value for 
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Table 3 

All-particle energy spectrum (d$s/dE) in units of (m 2 • sec • sr • GeV) -1 and cor- 
responding statistical (A stat ), total upper, (A+) and total lower (A_) errors and 
number of events (N ev ). The first four lines represent our data [TO] obtained for 
selection criteria R < 25m and 6 < 30°. 
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this significance of the bump depends somewhat on the energy range chosen to 
adjust the reference straight line in Fig. 9, but it lies in the range (3.5 — 4. 5) a. 

We conservatively included the systematic errors in this estimate, although 
they are not independent in the nearby points but correlated: the possible 
overestimation of the energy in one point cannot be followed by an underes- 
timation in the neighboring point if their energies are relatively close to each 
other. Systematic errors can change slightly the general slope of the spectrum 
but cannot imitate the fine structure and the existence of the bump. 

The results from Fig. 7 show that in the range of "bump" energy ( 7.4 • 10 7 
GeV) the systematic errors can not increase significantly the flux. To test this 
hypothesis more precisely we tested the reconstruction procedure for singular 
energy spectra with power indices 7 P = 1.5 and 7 P = 4.5 below and above the 
knee energy 7.4 • 10 7 GeV in the 5 • 10 7 — 10 s GeV energy range. Results are 
presented in Fig. 7 (star symbols) and show that there are no significant dis- 
crepancies of reconstructed spectra observed, which stems from high accuracy 
of energy reconstruction. 

The detected shower sample in the bump energy region did not reveal 
any discrepancies with the showers from adjacent energy bins within statis- 
tical errors neither with respect to reconstructed shower core coordinates, 
zenith angular and \ 2 distributions, nor with respect to £ = N ch /N^ distribu- 
tion. The only difference is that the average age of showers is increased from 
s = 0.88 ± 0.007 at E ~ 5 • 10 7 GeV up to s = 0.93 ± 0.01 at E ~ 10 8 GeV, 
instead of the monotonic shower age decrease with energy increment, which is 
observed at a low energy region (E = 3 • 10 6 — 5 • 10 7 GeV). 

It is necessary to note that some indications of the observed bump are also 
seen in KASCADE-Grande [6] (Fig. 9), TUNKA [4j and Tibet-Ill [27J data 
but with larger statistical uncertainties at the level of 1.5-2 standard devi- 
ations. Moreover, the locations of the bump in different experiments agree 
well with each other and with an expected knee energy for Fe-like primary 
nuclei according to the rigidity-dependent knee hypothesis [819"] . However, the 
observed width (~ 20% in energy) and height of the bump at the energy of 
(6 — 8) • 10 7 GeV, which exceeds by a factor of ~ 1.5 (~ 4 standard deviations) 
the best fit straight line fitting all points above 5 • 10 6 GeV in Fig. 9, are 
difficult to describe in the framework of the conventional model of cosmic ray 
origin [21]. 

As it will be shown below (Section 5, Fig. 10,11) the detected EAS charged 
particle (N c h) and muon size (A^) spectra [8f9] independently indicate the 
existence of this bump just for the obtained energies and as it follows from 
the behavior of shower age parameter versus shower size [8|9] , the bump at 
energy ~ 7.4 • 10 7 GeV is likely formed completely from Fe nuclei. 
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Fig. 10. EAS size spectra detected by the GAMMA facility (empty symbols) and 
corresponding expected spectra (filled symbols) computed in the framework of the 
SIBYLL2.1 interaction model and 2-component parametrization of primary spectra 
(11). The lines correspond to expected size spectra computed for each of primary 
nuclei. 



5 Possible origin of irregularities 



Irregularities of all-particle energy spectrum in the knee region are ob- 
served practically in all measurements [2|6|8] and are explained by both the 
rigidity-dependent knee hypothesis and contribution of pulsars in the Galac- 
tic cosmic ray flux [2Qf31f32] . Two these approaches approximately describe 
the all-particle spectrum in (1 — 100) • 10 6 GeV energy region. However, the 
observed bump in Fig. 9 at energies ~ 7.4 • 10 7 GeV both directly points out 
the presence of additional component in the primary nuclei flux and displays 
a very flat {p/ p ~ — 2) energy spectrum before a cut-off energy E c ~ 8 • 10 7 
GeV. 

It is known [8~|I9"] that rigidity-dependent primary energy spectra can not de- 
scribe quantitatively the phenomenon of ageing of EAS at energies (5 — 10) -10 7 
GeV that was observed in the most mountain-altitude experiments [8lll5|33j . 
It is reasonable to assume that an additional flux of heavy nuclei (Fe-like) 
is responsible for the bump at these energies. Besides, the sharpness of the 
bump (Fig. 9) points out at the local origin of this flux from compact objects 
(pulsars) [311132"] . 

The test of this hypothesis we carried out using parameterized inverse ap- 
proach [7|8]I9] on the basis of GAMMA facility EAS database and hypothesis 
of two-component origin of cosmic ray flux: 
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Fig. 11. The same as Fig. 10 for truncated EAS muon size spectra. 



F A (E) = <S>G(A)(E^(—y + P A (E)) (11) 
where Ph = Pne = Po = and 




The first term in the right hand side of expression (11) (so called Galactic 
component) is the power law energy spectra with rigidity-dependent knees at 
energies Ek = Er ■ Z and power indices 7 = 71 and 7 = 72 for E < E^ and 
E > Ek respectively, and the second term (so called "pulsar component") is an 
additional power law energy spectrum with cut-off energies E C;Fe and power 
indices 7 P = 71 )P and 7 P = 7 2p for E < E c Fe and E > E c Fe respectively. 

The scale factors &g(A) an d &p(A) along with particle rigidity E R , cut- 
off energy E C (A) and power indices 71, 72, 7i P , 72p were estimated using com- 
bined approximation method [7.8.9J for two examined shower spectra showed 
in Figs. 10,11 (empty symbols): EAS size spectra, dF/dN c h (Fig. 10) and EAS 
muon truncated size spectra, dF/dN^ (Fig. 11) detected by the GAMMA fa- 
cility using shower core selection criteria 9 < 30° and r < 50m [8,9j. We did 
not consider the H, He and O pulsar components to avoid a large number of 
unknown parameters and corresponding mutually compensative pseudo solu- 
tions [34] for the Galactic components. 

The folded (expected) shower spectra (filled symbols in Figs 10,11) were 
computed on the basis of parametrization (11) and CORSIKA EAS simulated 
data set [9115] for the A = H, He, O and Fe primary nuclei to evaluate the 
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Table 4 

Parameters of the primary energy spectra (11) derived from combined approxima- 
tions of detected shower spectra. The scale factors $>g,p{A) have units of (m 2 - s • 
sr • TeV)" 1 . The particle rigidity Er and cut-off energies E c are shown in PV and 
PeV units respectively. 



Par am. 


G-component 


P-component 




.102 ± .003 


— 


$(i?e) 


.094 ± .022 




HO) 


.032 ± .015 




$(Fe) 


.021 ± .006 


(.29 ± .08) • 10~ 7 


7i 


2.68 ± .005 


1.05 ± .5 


72 


3.29 ± .045 


4.5 ±.4 


Er 


2.59 ± .15 




E c ,Fe 




76.9 ± 1.5 



kernel functions of corresponding integral equations [HE]- The computation 
method, was completely the same as was performed in the combined approx- 
imation analysis [8119] . The initial values of spectral parameters for Galactic 
component were taken from [9f8] as well. In Figs. 10,11 we also presented the 
derived expected elemental shower spectra (lines) for primary H, He, O and 
Fe nuclei respectively. 

The parameters of two-component primary energy spectra (11) derived 
from the \ 2 goodness-of-fit test of shower spectra dF/dN ch and dF/dN^ are 
presented in Table 4. 

Resulting expected energy spectra Fa(E) for the Galactic H,He,0 and 
Fe nuclei (thin lines) along with the all-particle spectrum J2a^a(E) (bold 
line with shaded area) are presented in Fig. 12. The thick dash-dotted line 
corresponds to derived energy spectra of the additional Fe component (the 
second term in the right hand side of expression (11)). The all-particle energy 
spectrum obtained on the basis of the GAMMA EAS data and event-by-event 
multi-parametric energy evaluation method (Section 4, Fig.9) are shown in 
Fig. 12 (symbols) as well. 

It is seen, that the shape of 2-component all-particle spectrum (bold line 
with shaded area) calculated with parameters taken from the fit of EAS size 
spectra agrees within the errors with the results of event-by-event analysis 
(symbols) that points out at the consistency of applied spectral parametriza- 
tion (11) with GAMMA data. 

Notice that the flux of derived additional Fe component turned out to be 
about 0.5 — 0.6% of the total Fe flux for primary energies E > 10 6 GeV. This 
result agrees with expected flux of polar cap component [20J. 

The dependence of average nuclear mass number are presented in Fig. 13 
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Primary energy E (GeV) 



Fig. 12. All-particle primary energy spectrum (symbols) and expected energy spec- 
tra (lines and shaded area) derived from EAS inverse problem solution for p, He, O 
and Fe primary nuclei using 2-component parametrization (11). The thin lines are 
the energy spectra of Galactic H, He, O and Fe components. The thick dash-dotted 
line is an additional Fe component from compact objects. 




Primary energy E (GeV) 



Fig. 13. Average logarithm primary nuclei mass number derived from rigidity-de- 
pendent primary energy spectra [8 9j (dashed line) and 2-component model predic- 
tion (11) taking into account additional pulsar component (solid line). 
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for two primary nuclei flux composition models: one-component model, where 
the power law energy spectra of primary nuclei have rigidity-dependent knees 
at particle rigidity E R ~ 2500 GeV/Z [8"f9"] (so called Galactic component, 
dashed line) and two-component model (solid line), where additional pulsar 
(P) component was included according to parametrization (11) and data from 
Table 4 with very flat power index (ji p ~ 1) before cut-off energy E C} p e . The 
shaded area in Fig. 13 show the ranges of total (systematic and statistical) 
errors. 



6 Conclusion 



The multi-parametric event-by-event method (Sections 3,4) provides the 
high accuracy for the energy evaluation of primary cosmic ray nuclei cr(E) ~ 
10 — 15% regardless of the nuclei mass (biases< 5%) in the 3 — 200 PeV energy 
region. Using this method the all-particle energy spectrum in the knee region 
and above has been obtained (Fig. 9, Table 3) using the EAS database from 
the GAMMA facility. The results are obtained for the SIBYLL2.1 interaction 
model. 

The all-particle energy spectrum in the range of statistical and systematic 
errors agrees with the same spectra obtained using the EAS inverse approach 
[7116118] in the 3 — 200 PeV energy range. 

The high accuracy of energy evaluations and small statistical errors point 
out at the existence of an irregularity ('bump') in the 60 — 80 PeV primary 
energy region. 

The bump can be described by 2-component model (parametrization (11)), 
of primary cosmic ray origin, where additional (pulsar) Fe component are 
included with very flat power law energy spectrum (7^ ~ 1 ± 0.5) before 
cut-off energy E Ci F e (Fig. 13, Table 4). At the same time, the EAS inverse 
problem solutions for energy spectra of pulsar component forces the solutions 
for the slopes of Galactic component above the knee to be steeper (Table 4), 
that creates a problem of underestimation of all-particle energy spectrum in 
the range of HiRes g9] and Fly's Eye [28] data at E > 200 PeV. From this 
viewpoint the underestimation (10 — 15%) of the all-particle energy spectrum 
(bold solid line in Fig. 12) in the range of E > 200 PeV can be compensated 
by the expected extragalactic component [35J. 

Though we cannot reject the stochastic nature of the bump completely, our 
examination of the systematic uncertainties of the applied method lets us be- 
lieve that they cannot be responsible for the observed feature. The indications 
from other experiments mentioned in this paper provide the argument for the 
further study of this interesting energy region. 
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